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Abstract 

Point  clouds  are  one  of  the  most  primitive  and  fundamental  surface  representations.  A  popular  source  of  point 
clouds  are  three  dimensional  shape  acquisition  devices  such  as  laser  range  scanners.  Another  important  field 
where  point  clouds  are  found  is  in  the  representation  of  high- dimensional  manifolds  hy  samples.  With  the  increas¬ 
ing  popularity  and  very  broad  applications  of  this  source  of  data,  it  is  natural  and  important  to  work  directly 
with  this  representation,  without  having  to  go  to  the  intermediate  and  sometimes  impossible  and  distorting  steps 
of  surface  reconstruction.  A  geometric  framework  for  comparing  manifolds  given  by  point  clouds  is  presented 
in  this  paper.  The  underlying  theory  is  based  on  Gromov -Hausdorff  distances,  leading  to  isometry  invariant  and 
completely  geometric  comparisons.  This  theory  is  embedded  in  a  probabilistic  setting  as  derived  from  random 
sampling  of  manifolds,  and  then  combined  with  results  on  matrices  of  pairwise  geodesic  distances  to  lead  to  a 
computational  implementation  of  the  framework.  The  theoretical  and  computational  results  here  presented  are 
complemented  with  experiments  for  real  three  dimensional  shapes. 


1.  Introduction 

Point  clouds  are  one  of  the  most  primitive  and  fundamental 
manifold  representations.  One  of  the  most  popular  sources 
of  point  clouds  are  3D  shapes  acquisition  devices,  such  as 
laser  range  scanners,  with  applications  in  many  disciplines. 
These  scanners  provide  in  general  raw  data  in  the  form  of 
(noisy)  unorganized  point  clouds  representing  surface  sam¬ 
ples.  With  the  increasing  popularity  and  very  broad  appli¬ 
cations  of  this  source  of  data,  it  is  natural  and  important  to 
work  directly  with  this  representation,  without  having  to  go 
to  the  intermediate  step  of  fitting  a  surface  to  it  (step  that  can 
add  computational  complexity  and  introduce  errors).  See  for 
example  [4,  11,  13,  16,  21,  29,  30,  36,  37]  for  a  few  of  the 
recent  works  with  this  type  of  data.  Point  clouds  can  also  be 
used  as  primitives  for  visualization,  e.g.,  [5,  21,  40],  as  well 
as  for  editing  [44]. 

Another  important  field  where  point  clouds  are  found 
is  in  the  representation  of  high-dimensional  manifolds  by 
samples  (see  for  example  [2,24,41]).  This  type  of  high¬ 
dimensional  and  general  co-dimension  data  appears  in  al¬ 
most  all  disciplines,  from  computational  biology  to  image 
analysis  to  financial  data.  Due  to  the  extremely  high  dimen¬ 


sionality  in  this  case,  it  is  impossible  to  perform  manifold 
reconstruction,  and  the  task  needs  to  be  performed  directly 
on  the  raw  data,  meaning  the  point  cloud. 

The  importance  of  this  type  of  shape  representation  is 
leading  to  a  recent  increase  in  the  fundamental  study  of  point 
clouds  [1,  2,  9,  12,  17,  32,  33,  41]  (see  also  the  papers  men¬ 
tioned  in  the  first  paragraph  and  references  therein).  The 
goal  of  this  work,  ^  inspired  in  part  by  [14]  and  the  tools 
developed  in  [32,  41],  is  to  develop  a  theoretical  and  com¬ 
putational  framework  to  compare  shapes  (Riemannian  man¬ 
ifolds)  represented  as  point  clouds. 

As  we  have  mentioned,  a  variety  of  objects  can  be  repre¬ 
sented  as  point  clouds  in  Rd .  One  is  often  presented  with  the 
problem  of  deciding  whether  two  of  those  point  clouds,  and 
their  corresponding  underlying  objects  or  manifolds,  repre¬ 
sent  the  same  geometric  structure  or  not  (« object  recognition 
and  classification).  We  are  concerned  with  questions  about 


t  This  is  a  short  version  of  the  work  “A  Theoretical  and  Computa¬ 
tional  Framework  for  Isometry  Invariant  Recognition  of  Point  Cloud 
Data,”  by  Memoli  and  Sapiro,  submitted,  and  also  IMA  and  DTC 
Reports. 


© 


Memoli  and  Sapiro  /  Comparing  Point  Clouds 


the  underlying  unknown  structures  (objects),  which  need  to 
be  answered  based  on  discrete  and  finite  measures  taken  be¬ 
tween  their  respective  point  clouds.  In  greater  generality, 
we  wonder  what  is  the  structural  information  we  can  gather 
about  the  object  itself  by  exploring  the  point  cloud  which 
represents  it. 

Multidimensional  scaling  (MDS)  for  example  has  been 
used  to  approach  in  part  this  general  problem  of  object 
recognition.  Procedures  based  on  MDS  require  that  one 
first  computes  the  interpoint  distance  matrix  for  all  the 
members  of  the  point  cloud  (or  for  a  representative  selected 
sub- set  of  them).  If  one  is  interested  in  comparing  two 
different  objects,  the  problem  is  reduced  to  a  comparison 
between  the  corresponding  interpoint  distance  matrices.  If 
the  distance  we  use  is  the  Euclidean  one,  these  matrices 
only  provide  information  about  their  rigid  similarity,  and 
(assuming  they  have  the  same  size)  if  they  are  equal  (up 
to  a  permutations  of  the  indices  of  all  elements),  we  can 
only  conclude  that  there  exists  a  rigid  isometry  (rotation, 
reflection,  translation)  from  one  point  cloud  to  the  other. 
After  adding  compactness  considerations,  we  can  also  say 
something  about  the  true  underlying  (sampled)  objects. 
Being  a  bit  more  rigorous,  let  the  point  clouds  T?  C  Si  be 
^-coverings  of  the  surfaces  Si  in  R3,  for  i  =  1, 2  (this  will  be 
formally  defined  below).  Then  assuming  there  exists  a  rigid 
isometry  x  :  R3  ^  R3  such  that  x(05i)  =  T2,  we  can  bound 
the  Hausdorff  distance  (which  we  will  also  formally  define 
below)  between  x(Si)  and  S2  as  follows:  ^(-(x(5i ) .  1S2 )  < 
^r(A(‘S'i  ).t(Ti  ))  +  d^(tffP 1)^2)  +  ^r(<P2^2)  = 

dji(S i,Ti)  +  d‘K{rz(cy\),cS)2)  +  ^(^2,^2)  <  £1  +  0  +  82- 
And  of  course  the  same  kind  of  bound  holds  for  the  Haus¬ 
dorff  distance  between  the  points  clouds  once  we  assume 
the  underlying  continuous  objects  are  rigidly  isometric,  see 
§2.1  below. 

If  S\  and  S2  happen  to  be  isometric,  thereby  allowing 
for  bends  and  not  just  rigid  transformations,  we  wonder 
whether  we  will  be  able  to  detect  this  by  looking  at  (finite) 
point  clouds  T;  sampled  from  each  5).  This  problem  is  much 
harder  to  tackle.  We  approach  this  problem  through  a  proba¬ 
bilistic  model,  in  part  because  in  principle,  there  might  exist 
even  for  the  same  object,  two  different  samplings  that  look 
quite  dissimilar  (under  the  discrete  measures  we  can  cope 
with  computationally),  for  arbitrarily  fine  scales  (see  below). 

With  the  help  of  the  theory  presented  here  we  recast  these 
considerations  in  a  rigorous  framework  and  address  the  case 
where  the  distances  considered  to  characterize  each  point 
cloud  (object)  are  more  general.  We  concentrate  on  the  sit¬ 
uation  when  we  know  the  existence  of  an  intrinsic  notion 
of  distance  for  each  object  we  sample.  For  the  applications 
of  isometric  invariant  shape  (surfaces)  recognition,  one  must 
consider  the  distance  as  measured  by  paths  constrained  to 
travel  on  the  surface  of  the  objects,  better  referred  to  as 
geodesic  distance.  These  have  been  used  in  [14]  for  bending 
invariant  recognition  in  3D  (without  the  theoretical  founda¬ 


tions  here  introduced)  and  in  [17,  41]  to  detect  intrinsic  sur¬ 
face  dimensionality. 

In  this  paper,  the  fundamental  approach  used  for  isomet¬ 
ric  invariant  recognition  is  derived  then  from  the  Gromov- 
Hausdorff  distance  [19],  which  we  now  present.  If  two  sets 
(objects)  X  and  Y  are  subsets  of  a  common  bigger  metric 
space  (Z,  dz),  and  we  want  to  compare  X  to  Y  in  order  to  de¬ 
cide  whether  they  are/represent  the  same  object  or  not,  then 
an  idea  one  might  come  up  with  very  early  on  is  that  of  com¬ 
puting  the  Hausdorff  distance  between  them  (see  for  exam¬ 
ple  [10,  23]  for  an  extensive  use  of  this  for  shape  statistics 
and  image  comparison): 

dj{(X,Y)  =  max(sxipdz(x,Y),sxipdz(y,X)) 
xex  yer 

But,  what  happens  if  we  want  to  allow  for  certain  defor¬ 
mations  to  occur  and  still  decide  that  the  manifolds  are  the 
same?  More  precisely,  we  are  interested  in  being  able  to  find 
a  distance  between  metric  spaces  that  is  blind  to  isomet¬ 
ric  transformations  (“bends”).  This  will  permit  a  truly  ge¬ 
ometric  comparison  between  the  manifolds,  independently 
of  their  embedding  and  bending  position.  Following  [19], 
we  introduce  the  Gromov -Hausdorff  distance  between  Met¬ 
ric  Spaces 

dSx(X,Y)=  inf  4(Vn 

where  /  :  X  — >  Z  and  g  :  Y  — >  Z  are  isometric  embeddings 
(distance  preserving)  into  the  metric  space  Z.  It  turns  out  that 
this  measure  of  metric  proximity  between  metric  spaces  is 
well  suited  for  our  problem  at  hand  and  will  allow  us  to  give 
a  formal  framework  to  address  the  isometric  shape  recogni¬ 
tion  problem  (for  point  cloud  data).  However,  this  notion  of 
distance  between  metric  spaces  encodes  the  “metric”  dispar¬ 
ity  between  the  metric  spaces,  at  first  glance,  in  a  computa¬ 
tionally  impractical  way.  We  derive  below  new  results  that 
connect  this  notion  of  disparity  with  other  more  computa¬ 
tionally  appealing  expressions. 

Since  we  have  in  mind  specific  applications  and  scenar¬ 
ios  such  as  those  described  above,  and  in  particular  surfaces 
and  submanifolds  of  some  Euclidean  space  Rd ,  we  assume 
that  we  are  given  as  input  points  densely  sampled  from  the 
metric  space  (surface,  manifold).  This  will  manifest  itself  in 
many  places  in  the  theory  described  below.  We  will  present 
a  way  of  computing  a  discrete  approximation  (or  bound)  to 
)  based  on  the  metric  information  provided  by  these 
point  clouds.  Due  to  space  limitations,  the  proofs  are  omitted 
and  will  be  reported  elsewhere. 

2.  Theoretical  Foundations 

This  section  covers  the  fundamental  theory  behind  the  bend¬ 
ing  invariant  recognition  framework  we  develop.  We  use  ba¬ 
sic  concepts  of  metric  spaces,  see  for  example  [25]  for  a  de¬ 
tailed  treatment  of  this  and  [6,  19,  22,  26,  38,  39]  for  proofs 
of  Proposition  1  below. 
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Definition  1  (Metric  Space)  A  set  M  is  a  metric  space  if  for 
every  pair  of  points  x,  y  G  M  there  is  a  well  defined  function 
dM(x,y)  whose  values  are  non-negative  real  numbers,  such 
that  (a)  dM(x,y)  =  0^x  =  y,  and  (b)  du(x,y)  <  du(y,z)  + 
dM(z,x)  for  any  x,y,z  G  M.  We  call  dM'-MxM  — » R+  U {0} 
the  metric  or  distance.  For  clarity  we  will  specify  a  metric 
space  as  the  pair  (M^m)- 

Definition  2  (Covering)  For  a  point  x  in  the  metric  space 
(X,dx)  and  r  >  0,  we  will  denote  by  Bx(x,r)  the  set  {z  G 
X  :  dx(x,z)  <  r}.  For  a  subset  A  of  X ,  we  use  the  notation 
Bx(A,r)  =  Ufl£A5x(a,r).  We  say  that  a  set  C  C  X  is  an  R- 
covering  of  X  if  Bx(C,R)  =  X.  We  will  also  frequently  say 
that  the  set  A  is  a  ^-covering  of  X  if  A  constitutes,  for  some 
r  >  0,  a  covering  of  X  by  «-balls  with  centers  in  points  of  A. 

Definition  3  (Isometry)  We  say  the  metric  spaces  (X, dx) 
and  (Y,dy)  are  isometric  when  there  exists  a  bijective  map¬ 
ping  <|) :  X  — >  Y  such  that  dx(x  1,^2)  =  dy(<|>(jci),<|>(;c2))  for 
all  *1  ,*2  G  X.  Such  a  (|)  is  an  isometry  between  (X,dx)  and 

(Mr). 

Proposition  1 

1.  Let  (X.dx)i  (Y,dy)  and  (Z,dz)  be  metric  spaces  then 

d9x(X,r)  <  d9x(X,Z)  +  d9M(Z,Y). 

2.  If  (X ,  Y)  =  0  and  (X,dx),  (Y.  dy )  are  compact  metric 

spaces,  then  (X,dx)  and  (Y,dy)  are  isometric. 

3 .  Let  {x\ , . . . ,  xn  }  C  X  be  a  /^-covering  of  the  compact  met¬ 
ric  space  (X,dx)9  then  Jg^(X,  {*1, . . .  ,xn})  <  R. 

4.  For  compact  metric  spaces  (X,dx)  and  (Y,dy), 
i|D(X)-D(F)|  <d9K(X,Y)  <  imax(D(X),D(F)); 

where  D(X)  =  maxVJC/  ^xdx{x.x)  stands  for  the 
diameter  of  the  metric  space  X. 

5.  For  bounded  metric  spaces  (X,dx)  and  (Y,dy)  (x  G 
X,yGF), 

d$K  (X,  y )  =  inf  sup  \  | dx  (x,  \|/(y) )  -  dy  (y,  <|> (x))  \ 

x,y  Z 

V:  Y 

From  these  properties,  we  can  easily  prove  the  following 
important  result: 

Corollary  1  Let  X  and  Y  be  compact  metric  spaces.  Let 
moreover  Xm  be  a  r-covering  of  X  (consisting  of  m  points) 
and  Y m>  be  a  /  -covering  of  Y  (consisting  of  m  points).  Then 

\dw(X,Y)  -  ds^(Xm,Ymf)\ <  r  +  / 

We  can  then  say  that  if  we  could  compute  for 

discrete  metric  spaces  which  are  dense  enough  samplings 
of  the  continuous  underlying  ones,  that  number  would  be  a 
good  approximation  to  what  happens  between  the  continu¬ 
ous  spaces.  Currently,  there  is  no  computationally  efficient 
way  to  directly  compute  between  discrete  metric 

spaces  in  general.  This  forces  us  to  develop  a  roundabout 
path,  see  §2.2  ahead.  Before  going  into  the  general  case,  we 
discuss  next  the  application  of  our  framework  to  a  simpler 
but  important  case. 


2.1.  Intermezzo:  The  Rigid  Isometries  Case 

When  we  are  trying  to  compare  two  subsets  X  and  Y  of 
a  larger  metric  space  Z,  the  situation  is  less  complex.  The 
Gromov-Hausdorff  distance  boils  down  to  a  somewhat  sim¬ 
pler  Hausdorff  distance  between  the  sets.  In  more  detail, 

one  must  compute  d^^ld  (X ,  Y)  =  inf^>  d ^  (X ,  <|) (Y ) ) ,  where 
<1 > :  Z  ^  Z  ranges  over  all  self-isometries  of  Z.  If  we  know 
an  efficient  way  of  computing  inf<D^(  (X.d>(T)),  then  this 
particular  shape  recognition  problem  is  well  posed  for  Z, 
in  view  of  Corollary  1,  as  soon  as  we  can  give  guaran¬ 
tees  of  coverage.  This  can  be  done  in  the  case  of  sub¬ 
manifolds  of  Rd  by  imposing  a  probabilistic  model  on 
the  samplings  Xm  of  the  manifolds,  and  a  bound  on  the 
curvatures  of  the  family  of  manifolds.  In  more  detail  we 
can  show  that  P  {d1^  (X,  Xm)  >  8m^j  ~  as  m  |  00,  for 
/.  \l/k 

8m  >  (  )  ,  where  k  is  the  dimension  of  X,  see  Sec- 

tion  §3.2.  In  the  case  of  surfaces  in  Z  =  R  ,  <Y>  sweeps  all 
rigid  isometries,  and  there  exist  good  algorithms  which  can 
actually  solve  the  problem  approximately.  For  example,  in 
[18]  the  authors  report  an  algorithm  which  for  any  given 
0  <  a  <  1  can  find  <A>a  such  that 

43(Xm,$a(¥m,))  <  (8  +  a)  inWfc  (Xm,<J>(Ym<)) 

with  complexity  0(s4logs)  where  5  =  max(m,m/).  This 
computational  result,  together  with  our  theory,  makes  the 
problem  of  surface  recognition  (under  rigid  motions)  well 
posed  and  well  justified.  In  fact,  just  using  Corollary  1  and 
the  triangle  inequality,  we  obtain  a  bound  between  the  dis¬ 
tance  we  want  to  estimate  d^ngld  (X  ,F)  and  the  observed 
(computable)  value  d£(Xm,3>a(Y  m,)% 

\d%'rigid(X,Y)  — 43(Xm,$a(Ym-))|  < 

10  [r  +  r  +  d*^'rigid {X,Y)Y 

This  bound  gives  a  formal  justification  for  the  surface  recog¬ 
nition  problem.  To  the  best  of  our  knowledge,  this  is  the  first 
time  that  such  formality  is  shown  for  this  very  important 
problem,  both  in  the  particular  case  just  shown  and  in  the 
general  one  addressed  next. 


2.2.  The  General  Recognition  Case 

The  theory  introduced  by  Gromov  permits  to  address  the 
concept  of  (metric)  proximity  between  metric  spaces. 

When  dealing  with  discrete  metric  spaces,  as  those  aris¬ 
ing  from  samplings  or  coverings  of  continuous  ones,  it  is 
convenient  to  introduce  a  distance  between  them,  which  ul¬ 
timately  is  the  one  we  compute  for  point  clouds,  see  §3.4 
ahead.  For  discrete  metric  spaces  (both  of  cardinality  n ) 
(X  =  {*1,..., *„},dx)  and  (Y  =  {yi, . . .  ,y„},  dY)  we  define 
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the  distance: 

dj(X,Y)=  min  max  l-\dx{xi,Xj)-dY(yni,ynj)\  (1) 

TZE‘.Pn  1  <i.j<n  Z 

where  Yn  stands  for  the  set  of  all  n  x  n  permutations  of 
{ 1 , . . . ,  n}.  A  permutation  n  provides  the  correspondence  be¬ 
tween  the  points  in  the  sets,  and  \dx(xpXj)  —  dy ()%•, ) | 
gives  the  pairwise  distance/disparity  once  this  correspon¬ 
dence  has  been  assumed. 

It  is  evident  that  one  has  Y)  <  dy(K,  Y),by  virtue 

of  Property  5  from  Proposition  1 .  Moreover,  we  easily  derive 
the  following  easy  result,  whose  usefulness  will  be  made  ev¬ 
ident  in  §3. 

Corollary  2  Let  (X,dx)  and  (Y,dy)  be  compact  metric 
spaces.  Let  X  =  {x\ , . . . , xn}  C  X  and  Y  =  {y\ , . . .  ,yn}  C  Y, 
such  that  Bx(X,Rx)  =  X  and  By(Y,Ry)  =  Y  (the  point 
clouds  provide  Rx  and  Ry  coverings  respectively).  Then 

dSx (XJ)  <  RX  + Ry  +  d3  (X,  Y)  (2) 

with  the  understanding  that  dx  =  dx  \xxX  an<^  = 

dy  IyxY- 

Remark  1  This  result  tells  us  that  if  we  manage  to  find  cov¬ 
erings  of  X  and  Y  for  which  the  distance  dj  is  small,  then 
if  the  radii  defining  those  coverings  are  also  small,  the  un¬ 
derlying  manifolds  X  and  Y  sampled  by  these  point  clouds 
must  be  close  in  a  metric  sense.  Another  way  of  interpret¬ 
ing  this  is  that  we  will  never  see  a  small  value  of  Jj(X,Y) 
whenever  is  big,  a  simple  statement  with  prac¬ 

tical  value,  since  we  will  be  looking  at  values  of  dj,  which 
depend  on  the  point  clouds.  This  is  because,  in  contrast  with 
),  the  distance  dj  is  (approximately)  computable  from 
the  point  clouds,  see  §3.4. 

We  now  introduce  some  additional  notation  regarding  cov¬ 
erings  of  metric  spaces.  Given  a  metric  space  (X,dx),  the 

(r  s) 

discrete  subset  Nx  'n  denotes  a  set  of  points  {x\ , . . . ,  xn}  C  X 

such  that  (1)  Bx (NxlSJ ,  r)  =  X,  and  (2)  dx  { Xi,Xj )  >  5  when¬ 
ever  i  /  j.  In  other  words,  the  set  provides  a  coverage  and  the 
points  in  the  set  are  not  too  close  to  each  other  (the  coverage 
is  efficient). 

Remark  2  For  each  r  >  0  denote  by  N(r,X)  the  minimum 
number  of  closed  balls  of  radii  r  needed  to  cover  X.  Then, 
([38],  Chapter  10),  one  can  actually  show  that  the  class 
(M, of  all  compact  metric  spaces  X  whose  cover¬ 
ing  number  N(r,X)  are  bounded  for  all  (small)  positive  r 
by  a  function  N  :  (0,Ci)  — >  (0,  oo)  is  totally  bounded.  This 
means  that  given  p  >  0,  there  exist  a  finite  positive  inte¬ 
ger  fc(p)  and  compact  metric  spaces  X\ , . . .  ,Xk^  £  M  such 
that  for  any  X  £  M  one  can  find  i  £  {1, . . . , fc(p)}  such  that 
dcyj({X.Xi)  <  p.  This  is  very  interesting  from  the  point  of 
view  of  applications  since  it  allows  to  make  the  classifica¬ 
tion  problem  of  metric  spaces  in  a  well  justified  way.  For  ex¬ 
ample,  in  a  system  of  storage/retrieval  of  faces/information 
manifolds,  this  concept  permits  the  design  of  a  clustering 
procedure  for  the  objects. 


The  following  Proposition  will  also  be  fundamental  for  our 
computational  framework  in  §3,  leading  us  to  work  with 
point  clouds. 

Proposition  2  ([19])  Let  (X,dx)  and  (Y,dy)  be  any  pair  of 
given  compact  metric  spaces  and  let  r|  =  Jg^(X,T).  Also, 

(rs) 

let  NXn  =  {x\ , . . .  ,xn}  be  given.  Then,  given  a  >  0  there 
exist  points  {y^, . . .  ,y“}  C  Y  such  that 

!•  dj{N(xr'C{yi,  ■  ■  ■  ,y“})  <  (T|  +  a) 

2.  By  +  2(r|  +  a))  =Y 

3-  dY(yf,yJ)  >  s  -  2(r|  +  a)  for  i  ±  j. 

Remark  3  This  proposition  first  tells  us  that  if  the  metric 
spaces  happen  to  be  sufficiently  close  in  a  metric  sense,  then 
given  a  5- separated  covering  on  one  of  them,  one  can  find  a 
(s' -separated)  covering  in  the  other  metric  space  such  that  dj 
between  those  coverings  (point  clouds)  is  also  small.  This, 
in  conjunction  with  Remark  1,  proves  that  in  fact  our  goal 
of  trying  to  determine  the  metric  similarity  of  metric  spaces 
based  on  discrete  observations  of  them  is,  so  far,  a  (theoreti¬ 
cally)  well  posed  problem. 

Since  by  Tychonoff’s  Theorem  the  n- fold  product  space 

Y  x  ...  x  Y  is  compact,  if  5  —  2r|  >  c  >  0  for  some  con¬ 
stant  c,  by  passing  to  the  limit  along  the  subsequences  of 

{oc>0}  as  a  l  0  (if  needed)  above  one  can  as¬ 
sume  the  existence  of  a  set  of  different  points  {yi , . . .  ,yn}  C 

Y  such  that  dj({yi,. .  .,yn},N^'sJ)  <  ri,  min  i^jdriyijj)  > 
s-2r|  >  0,  and  By {{yi,--.,yn},r  +  2r\)  =  Y. 

Since  we  are  given  the  finite  sets  of  points  sampled  from 
each  metric  space,  the  existence  of  {yi , . . .  ,yn}  guaranteed 
by  Proposition  2  doesn’t  seem  to  make  our  life  a  lot  easier, 
those  points  could  very  well  not  be  contained  in  our  given  fi¬ 
nite  datasets  Xm  and  Ym/ .  The  simple  idea  of  using  a  triangle 
inequality  (with  metric  dj)  to  deal  with  this  does  not  work  in 
principle,  since  one  can  find,  for  the  same  underlying  space, 
two  nets  whose  dj  distance  is  not  small,  see  [7,  31].  Let  us 
explain  this  in  more  detail.  Assume  that  as  input  we  are  given 
two  finite  sets  of  points  Xm  and  Ym  on  two  metric  spaces, 
X  and  7,  respectively,  which  we  assume  to  be  isometric. 
Then  the  results  above  ensure  that  for  any  given  NXn  C  Xm 
there  exists  a C  Y  such  that  dj {Nx'sj , Ny'^ )  =  0.  How- 

(rs) 

ever,  it  is  clear  that  this  Ny'n  has  no  reason  to  be  con¬ 
tained  in  the  given  point  cloud  Yw.  The  obvious  idea  would 
be  try  to  rely  on  some  kind  of  independence  property  on 
the  sample  which  represents  a  given  metric  space,  namely 
that  for  any  two  different  covering  nets  N\  and  N2  (of  the 
same  cardinality  and  with  small  covering  radii)  of  X  the  dis¬ 
tance  dj  (N\ ,  N2)  is  also  small.  If  this  were  granted,  we  could 
proceed  as  follows:  d0(N^\N^])  <  + 

dj(N(py  ,N(yC)  =  0  +  small(r,  r,s,s),  where  small(r,  r,s, s) 
is  small  number  depending  only  on  r,  r ,  5  and  s.  The  property 
we  fancy  to  rely  upon  was  a  conjecture  proposed  by  Gromov 
in  [20]  (see  also  [42])  and  disproved  [7,31].  Their  coun- 
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terexamples  are  for  separated  nets  in  7Z1 .  It  is  not  known 
whether  we  can  construct  counterexamples  for  compact  met¬ 
ric  spaces,  or  if  there  exists  a  characterization  of  a  family  of 
^-points  separated  nets  of  a  given  compact  metric  space  such 
that  any  two  of  them  are  at  a  small  dj -distance  which  can  be 
somehow  controlled  with  n.  A  first  step  towards  this  is  the 
density  condition  introduced  in  [8]. 

If  counterexamples  do  not  exist  for  compact  metric 
spaces,  then  the  above  inequality  should  be  sufficient.  With¬ 
out  assuming  this,  we  give  below  an  argument  which  tackles 
the  problem  in  a  probabilistic  way.  In  other  words,  we  use 
a  probabilistic  approach  to  bound  dj  for  two  different  sam¬ 
ples  from  a  given  metric  space.  For  this,  we  pay  the  price  of 
assuming  the  existence  of  a  measure  which  comes  with  our 
metric  space.  On  the  other  hand,  probabilistic  frameworks 
are  natural  for  noisy  random  samples  of  manifolds  as  ob¬ 
tained  in  real  applications. 

2.3.  A  Probabilistic  Setting  for  Submanifolds  of  Rd 

We  now  limit  ourself  to  smooth  submanifolds  of  any  Rd, 
although  the  work  can  be  extended  to  more  general  metric 
spaces  (once  a  notion  of  uniform  probability  measure  is  in¬ 
troduced). 

Let  Z  be  a  smooth  and  compact  submanifold  of  Rd  with 
intrinsic  (geodesic)  distance  function  dz  We  can  now 
speak  more  freely  about  points  {zi}^L\  sampled  uniformly 
fromX.  For  any  measurable  C  C  Z,  P  (z.i  EC)  =  ||§|,  where 
a  (. B )  denotes  the  area  of  the  measurable  set  B  C  Z.  This  uni¬ 
form  distribution  can  be  replaced  by  other  distributions,  e.g., 
that  adapt  to  the  geometry  of  the  underlying  surface,  and  the 
framework  here  developed  can  be  extended  to  those  as  well. 

Let  Z  =  {z\ , . . . , Zn}  and  Zr  =  {z[ , . . . , zn}  be  two  discrete 
subsets  of  Z  (two  point  clouds).  For  any  permutation  n  G  ‘Yn 
and  i,je  {1, \dz(zi,Zj)  ~  dz{z%i,znj)\  <dz(zi,ZKi)  + 
dz(zj,zxj)  and  therefore  we  have 

df(Z,Z')  =  min  ma  xdz(zk,z%k)  >  dj(Z,Z')  (3) 

71GT„  k 

This  is  known  as  the  Bottleneck  Distance  between  Z  and  Z\ 
both  being  subsets  of  Z.  This  is  a  possible  way  of  measuring 
distance  between  two  different  samples  of  the  same  metric 
space. 

Instead  of  dealing  with  (3)  deterministically,  after  impos¬ 
ing  conditions  on  the  underlying  metric  space  Z,  we  derive 
probabilistic  bounds  for  the  left  hand  side.  We  also  make 
evident  that  by  suitable  choices  of  the  relations  among  the 
different  parameters,  this  probability  can  be  chosen  at  will. 
This  result  is  then  used  to  bound  the  distance  dj  between  two 
point  cloud  samples  of  a  given  metric  space,  thereby  lead¬ 
ing  to  the  bound  for  (a  quantity  related  to)  dj(Nx'sJ  ,Ny^) 
without  assuming  any  kind  of  proximity  of  the  nets  (and 
from  this,  the  bounds  on  the  original  Gromov-Hausdorff  dis¬ 
tance).  We  consider  Z  to  be  fixed,  and  we  assume  Z'  = 


{z[, . . .  ,zn}  to  be  chosen  from  a  set  Zm  C  Z  consisting  of 
m^>n  i.i.d.  points  sampled  uniformly  from  Z.  We  introduce 
the  Voronoi  diagram  V(Z)  on  Z,  determined  by  the  points 
in  Z  (see  for  example  [28]).  The  /- th  Voronoi  cell  of  the 
Voronoi  diagram  defined  by  Z  =  {z\ , . . .  ,Zn}  C  Z  is  given 

by  At  =  {zGZ:  dz(z,i,z)  <  min^-  dz(zj,z)}.  We  then  have 

Z  =  U5U4fc- 

We  first  want  to  find,  amongst  points  in  Zm,  n  different 
points  {zjj , . . . ,  Zin }  such  that  each  of  them  falls  inside  one 
Voronoi  cell,  {ztk  G  A^  for  k  =  1, . . .  ,n}.  We  provide  lower 
bounds  for  P(#(A^flZm)  >  1,  1  <k  <n),  the  probability 
of  this  happening.  We  can  see  the  event  as  if  we  collected 
points  inside  all  the  Voronoi  cells,  a  case  of  the  Coupon  Col¬ 
lecting  Problem ,  see  [15]:  we  buy  merchandise  at  a  coupon¬ 
giving  store  until  we  have  collected  all  possible  types  of 
coupons.  The  next  Lemma  presents  the  basic  results  we  need 
about  this  concept  [43]. 

Lemma  1  (Coupon  Collecting)  If  there  are  n  different 
coupons  one  wishes  to  collect,  such  that  the  probability  of 
seeing  the  k- th  coupon  is  p ^  (let  p  =  (p\ , . . . ,/?«)),  and  one 
obtains  samples  of  all  of  them  in  an  independent  way  then: 
The  probability  Pp(n,r)  of  having  collected  all  n  coupons 
after  r  trials  is  given  by 

/^(n,r)  =  l-5„^(-1)7^Wj  j  (4) 

where  the  symbol  Sn  means  that  we  consider  all  possible 
combinations  of  the  n  indices  in  the  expression  being  eval¬ 
uated.  (For  example  S3((pi  +  P2Y)  =  (pi  +  Pi)r  +  (pi  + 
PiY + {pi+ piY  ■) 

This  result  is  used  to  prove  the  following  fundamental 
probability  bounds: 

Theorem  1  Let  (Z,dz)  be  a  smooth  compact  submani¬ 
fold  of  Rd .  Given  a  covering  of  Z  and  a  number 

p  G  (0, 1),  there  exists  a  positive  integer  m  =  mn{p)  such 
that  if  Zm  =  {zk}™=i  is  a  sequence  of  i.i.d.  points  sam¬ 
pled  uniformly  from  Z,  with  probability  p  one  can  find 
a  set  of  n  different  indices  {i\  ,...,/„}  C  {1, . . .  ,m}  with 

d%(Nz*’{zi 

This  result  can  also  be  seen  the  other  way  around:  for 
a  given  m,  the  probability  of  finding  the  aforementioned 
subset  in  Zm  is  given  by  (4),  for  pz  defined  as  follows: 
plz  =  a(A*)/a(Z),  where  A/  is  the  /-th  Voronoi  cell  corre¬ 
sponding  to  N^\  1  <  /  <  n.  Moreover,  since  for  Zk  G  Nz  '^ 
Bz(zki^)  CAk  then  one  can  lower  bound  all  components  of 
pz.  In  practise  one  could  use  as  a  rule  of  thumb  m  ~  n\nn 
which  is  the  mean  waiting  time  (in  the  equiprobable  case) 
until  all  “coupons”  are  collected,  [15]. 

Corollary  3  Let  X  and  Y  compact  submanifolds  of  Rd.  Let 
Nx  „  be  a  covering  of  X  with  separation  5  such  that  for  some 
positive  constant  c,  s  —  2 Jg^(X,F)  >  c.  Then,  given  any 
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number  p  £  (0, 1),  there  exists  a  positive  integer  m  =  mn(p) 
such  that  if  Ym/  =  {yk}k=i  is  a  sequence  of  i.i.d.  points  sam¬ 
pled  uniformly  from  Y,  we  can  find,  with  probability  at  least 
p ,  a  set  of  n  different  indices  {4 , . . . ,  4}  C  {1, . . . ,  m}  such 
that  dj(N^,{yh , . . .  <  3  dm(XJ)  +  r. 

This  concludes  the  main  theoretical  foundation  of  our  pro¬ 
posed  framework.  We  have  shown  that  dj  is  a  good  approxi¬ 
mation  of  the  Gromov-Hausdorff  distance  between  the  point 
clouds,  in  a  probabilistic  sense.  Now,  we  must  devise  a  com¬ 
putational  procedure  which  allows  us  to  actually  find  the 
subset  {y/j , . . . ,  ytn }  inside  the  given  point  cloud  Ym  when 
it  exists,  or  at  least  find  it  with  a  large  probability.  Note  that 
in  practise  we  can  only  access  metric  information,  that  is, 
interpoint  distances.  Point  positions  cannot  be  assumed  to 
be  accessible  since  that  would  imply  knowing  the  (isome¬ 
try)  transformation  that  maps  X  into  Y.  A  stronger  result 
should  take  into  account  possible  self-isometries  of  X  (F), 
which  would  increase  the  probability  of  finding  a  net  which 
achieves  small  dj  distance  to  the  fixed  one.  Next  we  present 
such  a  computational  framework. 

3.  Computational  Foundations 

There  are  a  number  of  issues  that  must  be  addressed  in  or¬ 
der  to  develop  an  algorithmic  procedure  from  the  theoretical 
results  previously  presented.  These  are  now  addressed. 

3.1.  Initial  Considerations 

In  practise  we  have  as  input  two  independent  point  clouds 
Xm  and  Ym/  each  of  them  composed  of  i.i.d.  points  sam¬ 
pled  uniformly  from  X  and  F,  respectively.  We  fix  a  number 
n  <  min  (m,  m)  and  construct  good  coverings  Nx'„  of  X  and 

NynS  ^  of  F.  Actually,  r,s,r'  and  s'  all  depend  on  n ,  and  we 
should  choose  n  such  that  r  and  r'  are  small  enough  to  make 
our  bounds  useful,  see  the  additional  computations  below. 
Details  on  how  we  construct  these  coverings  are  provided  in 
Section  §3.3. 

It  is  convenient  to  introduce  the  following  additional  no¬ 
tation.  For  q  ON,  let  {1  :  q}  denote  the  set  {1  also 

for  a  set  of  points  Zq  =  {zk }|=1  and  for  a  set  of  1  <  u  <  q 
indices  4  =  {h ,  •  •  • ,  4}  C  {1  :  q},  let  Zq[Iu\  denote  the  set 

\z,ix , .  .  . 

Corollary  3  suggests  that  in  practise  we  compute  the  sym¬ 
metric  expression 

%  (Xm,Ym)  =  (5) 

max  (  min  dj(Np\V„[J„]),  min  dj(Nyns \X44i])) 
\Jn  C  { l:m}  ’  In  C  f  t:m}  ’  J 

which  depends  not  only  on  Xm  and  Ym>  but  also  on  specified 

(  V  5^  ( s ?  1 

covering  nets  NXn  and  NYn  .  However  we  prefer  to  omit 
the  dependence  in  the  list  of  arguments  to  keep  the  notation 
simpler. 


Then,  we  know  that  with  probability  at  least  Ppx  (n,m)  x 
PpY(n,m  )  we  have  (we  assume  Xm  to  be  independent 
from  Ym/)  Jgr(Xm,  Ym/)  <  3 dg^(X,Y)  +max(r,r/).  More¬ 
over,  in  some  precise  sense  Ym/)  upper  bounds 

Jgy^(Xm,  Ym/),  something  we  need  to  require  otherwise  we 
would  have  solved  one  problem  to  gain  another,  and  that  im¬ 
plies  (Corollary  1)  a  similar  upper  bound  for  dc;^(X,F). 

In  fact,  for  any  4  C  {1  :  m} 

dso-c0^m,Ym /)  <  ^giH;(Xm, Xm[4])  +^giH;(Xm[4],  Ym/) 

<  d9M(Xm,Xm[In})  +d3M(Xm[In]AY^ 

+  dSx(N^s%Ym,) 

<  d^(Xm,X4ln})  +  d3{X4ln},N^'/*)  +  r' 

(r1  s1 ) 

Now,  considering  4  such  that  dj  (X  m  [In],NYn  )  = 

min/„C{l:m}^(Yn)’Xm[/"])'  We  find  dSJt(Xm,Ym/)  < 
dj^(Xm,Xm[In\) dgr(Xm,Y mi) r  . 

Symmetrically,  we  also  obtain  for  Jn  such  that 

dj(Ym[Jn],N(x»)  =  miny„c{l:m'}  dj(N{xry,Ym,[Jn  ]) 

d3 jc(Xm,  Ym<)  <  Ymf[Jn])  +d3r(Xm, Ymf)  +r 

Hence,  Jg^(Xm,  Ym/)  <  %(Xm,  Ym/)  + 

min  (Xm,  Xw[4] ) ,  ( Yw/ ,  Ym>  [4] ))  +  max(r,  r) . 

Let  Ax  =  Xm[4|J  and  A y  =  d^(Ym/,Ym/[Jn]). 

The  computational  procedure  we  infer  is:  If  4y(Xm,  Ym/) 
is  “large”we  then  know  that  Jgy^(X,T)  must  also  he 
“large” with  high  probability.  On  the  other  hand,  if 
%(Xm,Ym/)  is  “small”  and  min  (Ax,  Ay)  is  also  “small” 
then  dc;yt(X,r)  must  also  be  “small.” 

3.2.  Working  with  Point  Clouds 

First,  all  we  have  is  a  finite  sets  of  points,  point  clouds,  sam¬ 
pled  from  each  metric  space,  and  all  our  computations  must 
be  based  on  these  observations  alone.  Since  we  made  the  as¬ 
sumption  of  randomness  in  the  sampling  (and  it  also  makes 
sense  in  general  to  model  the  problem  in  this  way,  given  the 
way  shapes  are  acquired  by  a  scanner  for  example),  we  must 
relate  the  number  of  acquired  points  to  the  coverage  prop¬ 
erties  we  wish  to  have.  In  other  words,  and  following  our 
theory  above,  we  would  like  to  say  that  given  a  desired  prob¬ 
ability  p  and  a  radius  r,  there  exists  a  finite  m  such  that  the 
probability  of  covering  all  the  metric  space  with  m  balls  (in¬ 
trinsic  or  not)  of  radius  r  centered  at  those  m  random  points 
is  at  least  p.  This  kind  of  characterizations  are  easy  to  deal 
with  in  the  case  of  submanifolds  of  Rd ,  where  the  tuning 
comes  from  the  curvature  bounds  available.  For  this  we  fol¬ 
low  arguments  from  [32].  Let  Z  be  a  smooth  and  compact 
submanifold  of  Rd  of  dimension  k.  Let  Zm  C  Z  consist  of 
m  i.i.d.  points  uniformly  sampled  from  Z.  Let  K  be  an  upper 
bound  for  the  sectional  curvatures  of  Z.  Then  we  can  prove 
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that  for  a  sequence  rm  — »  0  such  that  rm  >  ^  for  large  m, 

P  (4  (Z,ZW)  >  ~ 

Then,  since  one  can  also  prove,  [32],  that  for  any  z  G  Z, 
5  >  0  small,  B(z,8)  (TZ  C  Bz(z,Ck8),  for  some  constant 
Cf  >  1  depending  only  on  metric  properties  of  Z  (curvatures 
and  diameter),  we  also  find  P  Zm)  >  rm^J  ~ 

This  relation  gives  us  some  guidance  regarding  how  many 
points  we  must  sample  in  order  to  have  a  certain  covering  ra¬ 
dius,  or  to  estimate  the  covering  radius  in  terms  of  m.  More 
precise  estimates  can  be  found  in  the  reference  mentioned 
above.  The  important  point  to  remark  is  that  this  kind  of  re¬ 
lations  should  hold  for  the  family  of  shapes  we  want  to  work 
with,  therefore,  once  given  bounds  on  the  curvatures  and  di¬ 
ameters  which  characterize  the  family,  one  can  determine  a 
precise  probabilistic  covering  relation  for  it.  We  leave  the 
exploitation  of  this  idea  for  future  work. 

Given  the  natural  number  n<^m( or  eventually  5  >  0),  we 
use  the  procedure  described  in  §3.3  below  to  find  ^-points 
from  Zm  which  constitute  a  covering  of  Zm  of  the  given  car¬ 
dinality  n  (or  of  the  given  separation  s)  and  of  a  resulting 

irs) 

radius  r.  We  denote  this  set  by  Nz  n  C  Zm. 

3.3.  Finding  Coverings 

In  order  to  find  the  coverings,  we  use  the  well  known  Far¬ 
thest  Point  Sampling  (FPS)  procedure,  which  we  describe 
next.  Suppose  we  have  a  dense  sampling  Zm  of  the  smooth 
and  compact  submanifold  of  Rd  ( Z,dz )  as  interpreted  by  the 
discussion  above.  We  want  to  simplify  our  sampling  and  ob¬ 
tain  a  well  separated  covering  net  of  the  space.  We  also  want 
to  estimate  the  covering  radius  and  separation  of  our  net.  It  is 
important  to  obtain  subsets  which  retain  as  best  as  possible 
the  metric  information  contained  in  the  initial  point  cloud 
in  order  to  make  computational  tasks  more  treatable  without 
sacrificing  precision. 

We  first  show  a  procedure  to  sample  the  whole  space  Z. 
Fix  n  the  number  of  points  we  want  to  have  in  our  simplified 
point  cloud  3V  We  build  yn  recursively.  Given  T„_i,  we  se¬ 
lect  p  G  Z  such  that  dz(p^n)  =  maxzGz dz(z,  Tw_i)  (here 
we  consider  of  course,  geodesic  distances).  There  might  ex¬ 
ist  more  than  one  point  which  achieves  the  maximum,  we 
either  consider  all  of  them  or  randomly  select  one  and  add  it 
to  OVi.  This  subsampling  procedure  has  been  studied  and 
efficiently  implemented  in  [34]  for  the  case  of  surfaces  rep¬ 
resented  as  point  clouds. 

Let  us  now  assume  that  the  discrete  metric  space  (Zm,  dz) 
is  a  good  random  sampling  of  the  underlying  (Z,dz)  in  the 
sense  that  d^(Z,  Zm)  <  r  with  probability  pr,m,  as  discussed 
in  Section  §3.2.  We  then  want  to  simplify  Zm  in  order  to 
obtain  a  set  with  n  points  which  is  both  a  good  subsam¬ 
pling  and  a  well  separated  net  of  Z.  We  want  to  use  our  n 
sampled  points  in  the  best  possible  way.  We  are  then  led  to 


using  the  construction  discussed  above.  Choose  randomly 
one  point  p\  G  Zm  and  consider  Tj  =  {p\}.  Run  the  proce¬ 
dure  FPS  until  n  —  1  other  points  have  been  added  to  the  set 
of  points.  Compute  now  rn  =  m2LXq^mdz(q,  T«).  Then,  also 
with  probability  pr?m,  is  a  (r  +  rn) -covering  net  of  Z  with 
separation  sn,  the  resulting  separation  of  the  net.  Following 
this,  we  now  use  the  notation  N^n  . 

We  use  a  graph  based  distance  computation  following  [3], 
or  the  exact  distance,  which  can  be  computed  only  for  certain 
examples  (spheres,  planes).  We  could  also  use  the  techniques 
developed  for  triangular  surfaces  in  [27],  or,  being  this  the 
optimal  candidate,  the  work  on  geodesics  on  (noisy)  point 
clouds  developed  in  [32]. 

3.4.  Additional  Implementation  Details 

In  this  section  we  conclude  the  details  on  the  implementa¬ 
tion  of  the  framework  here  proposed.  The  first  step  of  the 
implementation  is  the  computation  of  dj  and  subsequently 
c/gr,  which  from  the  theory  we  described  before,  bounds  the 
Gromov-Hausdorff  distance. 

We  have  implemented  a  simple  algorithm.  Considering 
the  matrix  of  pairwise  geodesic  distances  between  points  of 
Xm,  we  need  to  determine  whether  there  exists  a  submatrix 
of  the  whole  distance  matrix  corresponding  to  Xm  which  has 

a  small  dj  distance  to  the  corresponding  matrix  of  a  given 

*y  s>\ 

Nyn  .  We  select  this  latter  net  as  the  result  of  applying  the 
FPS  procedure  to  obtain  a  subsample  consisting  of  n  points, 
where  the  first  two  points  are  selected  to  be  at  maximal  dis¬ 
tance  from  each  other.  To  fix  notation,  let  Xm  =  {x\ , . . .  ,xm} 

and  Nyn's  ^  =  {yj] , . . .  ,y7n}.  We  then  use  the  following  al¬ 
gorithm. 

(k  =  1,2)  Choose  Xi]  and  v;2  such  that  \dx(xix  ,jc/2)  — 
dy  (; yp ,  yj2 )  |  is  minimized. 

(k  >  2)  Let  xik  ,  G  Xm  be  such  that  (x4+1)  = 

mini<j,<m ek+i  ixii )  where  eM  (xu)  = 

maxi<r<*  \dx(xinXir)-dY{yjk+l  ,yjr)\. 

We  stop  when  n  points,  {xi] ,  jc,-2  , . . .  }  have  been  selected, 
and  therefore  a  distance  submatrix  ((dx(xiu,Xiv)))”v= i,  is 

obtained.  Since  we  can  write  dj  ( {xi: , . . . ,  } ,  Nyn's  ^)  = 

\ maxi <k<n maxj <t<k-\\dx (xik ,xit )  -  dY(yjk,yj,)\  = 
2 maxi <k<„ ekixir)  we  then  see  that  with  our  algorithm  we 
are  minimizing  the  error  point- wise. 

Of  course,  we  now  use  the  same  algorithm  to  compute  the 
other  half  of  dj.  This  algorithm  is  not  computationally  opti¬ 
mal.  We  are  currently  studying  computational  improvements 
along  with  error  bounds  for  the  results  provided  by  the  algo¬ 
rithms. 

4.  Examples 

We  now  present  experiments  that  confirm  the  validity  of  the 
theoretical  and  computational  framework  introduced  in  pre- 
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vious  sections.  In  the  future,  we  plan  to  make  these  experi¬ 
ments  more  rigorous,  including  concepts  of  hypothesis  test¬ 
ing.  As  a  simplification,  for  our  experiments  we  have  only 
computed  d ^  neglecting  the  other  terms  (see  §3.1)  which 
would  provide  a  estimative  of  the  Gromov-Hausdorff  prox¬ 
imity  between  the  shapes. 

We  complemented  the  more  complex  data  (as  presented 
below)  with  simple  shapes:  (1)  A  plane,  Pn  =  [—  ^=,  ~^]2 

and  (2)  A  sphere,  S  =  {x  £  Rd  :  ||x||  =  1}. 

We  first  test  our  framework  when  X  and  Y  are  isometric. 
We  first  consider  X  =  Y  and  see  whether  we  make  the  right 
decision  based  on  the  discrete  (random)  measurements.  Let 
Xm  and  Ym  be  two  independent  sets  composed  of  m  indepen¬ 
dent,  uniformly  distributed  random  points  on  X.  In  the  case 
of  the  sphere  we  generated  this  uniformly  distributed  sample 
points  using  the  method  of  Muller,  see  [35].  We  consider  X 
to  be  either  the  plane  Pn  or  the  sphere  S  as  defined  above. 
Given  n ,  from  Xm  and  Y m,  and  using  the  FPS  procedure,  we 
construct  NymJl  and  Ny m,n  (we  omit  the  supraindices  since 
we  won’t  use  the  values  of  covering  radius  and  separation), 
and  look  for  a  metric  match  inside  Xm  and  Y m,  respectively, 
following  the  algorithm  described  in  §3.4  for  the  computa¬ 
tion  of  dgr(Xm,  Ym).  (Recall  that  actually  d^(Xm,  Ym)  de¬ 
pends  on  n ,  see  its  definition  (5).)  For  each  dataset  we  tested 
for  values  of  m  £  M  =  {500, 600, . . . ,  2000}  and  n  E  N  = 
{5, 10, 15, ... ,  100},  and  obtained  the  results  reported  below. 
In  Table  1  we  show  the  values  of  dj  for  selected  values  of 
m  and  n.  As  expected,  the  values  of  d<$  are  small  compared 
to  D  (Pn)  =  D  (S)  =  71  (see  below  for  the  corresponding  val¬ 
ues  when  comparing  non-isometric  shapes).  In  Figure  1  (first 
two  figures)  we  show  a  pseudocolor  representation  of  the  re¬ 
sults  for  dj. 

We  now  proceed  to  compare  shapes  that  are  not  isomet¬ 
ric,  starting  with  X  =  P%  (a  plane)  and  Y  =  S  (a  sphere).  In 
this  case  we  expect  to  be  able  to  detect,  based  on  the  finite 
point  clouds,  that  d ^  is  large.  Table  1  (see  also  last  two  fig¬ 
ures  of  Figure  1),  shows  the  results  of  a  simulation  in  which 
we  compared  the  sphere  S  and  the  plane  Pn,  varying  the  net 
sizes  and  the  total  number  of  points  uniformly  sampled  from 
them.  The  experiments  have  been  repeated  100  times  to  pro¬ 
duce  this  table,  and  the  reported  values  consist  of  the  mean 
of  these  100  tests,  as  well  as  their  maximum  (the  correspond¬ 
ing  deviation  was  1.72  x  10-2).  As  expected,  the  values  are 
larger  than  when  comparing  plane  against  plane  or  sphere 
against  sphere. 

We  conclude  the  experiments  with  real  (more  complex) 
data.  We  have  4  sets  of  shapes  (the  datasets  were  kindly 
provided  to  us  by  Prof.  Kimmel  and  his  group  at  the  Tech- 
nion),  each  one  with  their  corresponding  bends.  We  ran  the 
algorithm  N  =  6  times  with  n  =  70,  m  =  2000,  using  the 
4  nearest  neighbors  to  compute  the  geodesic  distance  using 
the  isomap  engine.  The  data  description  and  results  are  re¬ 
ported  in  Table  2.  We  note  not  only  that  the  technique  is  able 
to  discriminate  between  different  object,  but  as  expected,  it 


n\m 

500 

900 

1500 

1900 

5 

0.036793 

0.015786 

0.018160 

0.0074027 

25 

0.041845 

0.050095 

0.026821 

0.031019 

45 

0.081975 

0.042198 

0.038990 

0.036376 

65 

0.068935 

0.052482 

0.035718 

0.031512 

85 

0.077863 

0.038660 

0.036009 

0.036894 

n\m 

500 

900 

1500 

1900 

5 

0.013282 

0.013855 

0.010935 

0.013558 

25 

0.082785 

0.043617 

0.033095 

0.033592 

45 

0.074482 

0.067096 

0.057161 

0.040727 

65 

0.079456 

0.076762 

0.049503 

0.043405 

85 

0.083577 

0.083344 

0.058094 

0.054144 

n\m  500  1000  1500  2000 


10 

1.839 

X 
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1.902 

X 

10-1 

1.931 

X 

10-1 

1.942 

X 

10“ 

25 

1.834 

X 

10_1 

1.908 

X 

10_1 

1.920 

X 

10_1 

1.944 

X 

10“ 

50 

1.818 

X 

10_1 

1.899 

X 

10_1 

1.925 

X 

10_1 

1.933 

X 

10“ 

75 

1.873 

X 

10-1 

1.882 

X 

10-1 

1.936 

X 

10-1 

1.939 

X 

10“ 

100 

1.846 

X 

10_1 

1.913 

X 

10_1 

1.924 

X 

10_1 

1.936 

X 

10“ 

Table  1:  Table  with  values  of  dj  for  a  plane  (top),  a  sphere 
(middle),  and  a  plane  against  a  sphere  (bottom). 


doesn’t  get  confused  by  bends.  Moreover,  the  distances  be¬ 
tween  a  given  object  and  the  possible  bends  of  another  one 
are  very  similar,  as  it  should  be  for  isometric  invariant  recog¬ 
nition. 


5.  Conclusions 

A  theoretical  and  computational  framework  for  comparing 
manifolds  (metric  spaces)  given  by  point  clouds  was  intro¬ 
duced  in  this  paper.  The  theoretical  component  is  based  on 
the  Gromov-Hausdorff  distance,  which  has  been  extended 
and  embedded  in  a  probabilistic  framework  to  deal  with 
point  clouds  and  computable  distances.  Examples  support¬ 
ing  this  theory  were  provided. 

We  are  currently  working  on  improving  the  computational 
efficiency  of  the  algorithm,  performing  additional  experi¬ 
ments,  and  in  particular,  comparing  high  dimensional  point 
clouds  with  data  from  image  sciences  and  neuroscience.  This 
as  well  as  the  proofs  of  the  theorems  in  this  paper  will  be 
published  elsewhere. 
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Figure  1:  First  two  figures:  Graphic  visualization  of  the  results  for  the  plane  Pn  (on  the  left),  and  the  sphere  S  (on  the  right).  Red  corresponds 
to  low  values  of  dj  and  blue  larger  values.  On  the  horizontal  axis,  from  left  to  right  we  have  increasing  values  of  m,  while  on  the  vertical  axis, 
n  increases  going  upwards.  Observe  how  the  distance  increases  for  fixed  m  as  n  increases  in  accordance  with  the  fact  that  we  have  less  freedom 
to  choose  the  n  points  from  the  given  m.  Third  figure:  Graphic  visualization  of  the  results  for  the  comparison  between  the  plane  Pn  and  the 
Sphere  S.  Red  corresponds  to  low  values  of  dj  and  blue  larger  values.  On  the  horizontal  axis,  from  left  to  right  we  find  increasing  values  of  m, 
and  on  the  vertical  axis,  n  increases  going  upwards.  Fourth  figure:  Plot  of  the  values  of  d<j  obtained  against  n,  the  size  of  the  FPS  net,  with 
m  =  2000.  (This  is  a  color  figure.) 
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5.95  X  10“ 1  5.95  X  10“ 1  3.482  3.464  1.12  x  10“ 


T 


4.19  Xl0_1  4.19  Xl0_1  3.31  3.29  1.62  xlO-1  1.59  xlO-1 


4.25  Xl0_1  4.25  Xl0_1  3.31  3.29  1.56  X  10_1  1.15xlO_i  5.53x10 


1  i  i  s  v  in- 1 


4.16  Xl0_1  4.16  Xl0_1  3.30  3.28  1.65  xlO-1  1.62  xlO-1  4.85  X  10-2  5.77  X  10-2 


Diameters  1.223  1.223  6.996  6.960  6.1  X  10-2  6.8  X  10-2  3.86  X  10_1  3.73  xlO-1  3.91x10 


r>—  1  n?  v  1  1 


Figure  2:  Comparison  results  for  the  complex  objects.  The  number  of  points  per  model  are  indicated  in  the  first  row  under  the 
corresponding  figure. 
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